###########################################################
# Normalization of Censorship: Evidence from China
# Replication Files (Experimental Studies, Main Paper)
#
# Author: Tony Zirui Yang
# Date: January 19, 2024
#
# The lines below reproduce the main Figures and results
# related to the experimental studies reported in the paper
###########################################################

rm(list=ls())

# Packages ----
library("dplyr")
library("ggplot2")
library("tidyverse")
library("rstatix")
library("ggpubr")
library("gridExtra")
library("mediation")
library("data.table")

# Load the data
Experiment <- fread("Experiment.csv")
# Study 1 data
Study1 <- Experiment[which(Experiment$Study == 1),]
# Study 2 data
Study2 <- Experiment[which(Experiment$Study == 2),]


## ------------------------------
## Figure 3 
## ------------------------------

# Outcome variables
formula <- c("Censor_Support ~ ", "Regime_Overall ~ ", "Regime_Central ~ ",
             "Regime_Local ~ ", "Regime_Protest ~ ")

# Study 1 Results (Difference in Means)
lm_list_s1 <- lapply(formula, function(f) lm(paste0(f, "Treatment"), data = Study1))
names(lm_list_s1) <- c("lm_censor_s1", "lm_overall_s1", "lm_central_s1", "lm_local_s1", "lm_protest_s1")

# Study 2 Results (Difference in Means)
lm_list_s2 <- lapply(formula, function(f) lm(paste0(f, "Treatment"), data = Study2))
names(lm_list_s2) <- c("lm_censor_s2", "lm_overall_s2", "lm_central_s2", "lm_local_s2", "lm_protest_s2")

# Two Studies Combined Results (Difference in Means)
lm_list_c <- lapply(formula, function(f) lm(paste0(f, "Treatment"), data = Experiment))
names(lm_list_c) <- c("lm_censor_c", "lm_overall_c", "lm_central_c", "lm_local_c", "lm_protest_c")

# Making Figure 3
# Variable Names
Terms <- factor(rep(c("Support for\nCensorship Apparatus", "Regime Support:\nOverall Satisfaction",
                      "Regime Support:\nCentral Government", "Regime Support:\nLocal Government",
                      "Willingness to Protest"),each = 3),
                levels=rev(c("Support for\nCensorship Apparatus", "Regime Support:\nOverall Satisfaction",
                             "Regime Support:\nCentral Government", "Regime Support:\nLocal Government",
                             "Willingness to Protest")), ordered=TRUE)
# Combining Models into a list
models <- c(lm_list_s1, lm_list_s2, lm_list_c)
models <- models[c("lm_censor_s1", "lm_censor_s2", "lm_censor_c",
                   "lm_overall_s1", "lm_overall_s2", "lm_overall_c",
                   "lm_central_s1", "lm_central_s2", "lm_central_c",
                   "lm_local_s1", "lm_local_s2", "lm_local_c",
                   "lm_protest_s1", "lm_protest_s2", "lm_protest_c")]
# Point Estimates
Estimates <- sapply(models, function(x) coef(summary(x))["Treatment","Estimate"])
# Standard Errors
Std_Error <- sapply(models, function(x) coef(summary(x))["Treatment","Std. Error"])
# 95% Confidence Intervals
CI95 <- t(sapply(models, function(x) confint(x, level = 0.95)["Treatment",]))
# 90% Confidence Intervals
CI90 <- t(sapply(models, function(x) confint(x, level = 0.90)["Treatment",]))
# Legend
Study <- rep(c("Study 1", "Study 2", "Two Studies Combined"), 5)
# Create a data frame for the Figure 3
Figure.3 <- data.frame(Terms, Estimates, Std_Error, Study, CI95, CI90)
# Rename data frame column names
Figure.3 <- Figure.3 %>% 
  rename("conf.low_95" = "X2.5..",
         "conf.high_95" = "X97.5..",
         "conf.low_90" = "X5..",
         "conf.high_90" = "X95..")
# Adjust positions between bars
pd = position_dodge2(width = 0.4, reverse = TRUE)
# ggplot
ggplot(Figure.3, aes(x=Terms, y = Estimates, color= Study, shape = Study)) + 
  # Line at 0 intercept
  geom_hline(yintercept = 0, colour = gray(1/2), lty = 2) + 
  # Point estimates
  geom_point(position = pd, size = 2) + 
  # Confidence intervals
  geom_linerange(aes(ymin = conf.low_90, ymax = conf.high_90), lwd = 0.9, position = pd) + 
  geom_linerange(aes(ymin = conf.low_95, ymax = conf.high_95), lwd = 0.5, position = pd) + 
  scale_y_continuous(limits = c(-0.5, 0.5)) +
  # Labels
  labs(y="Treatment Effects", x="Outcome Variables") +  
  # Rotate the plot
  coord_flip() +  
  # Nicer theme 
  theme_bw() + 
  theme(axis.text = element_text(size = 12, face = "bold"),
        axis.title = element_text(size=12),
        axis.text.y = element_text(size = 12)) 

# Saving the ggplot
ggsave(file = "Figures_Paper/Figure_03.pdf",   # The directory you want to save the file in
       width = 9, # The width of the plot in inches
       height = 5) # The height of the plot in inches


# Adjusting p-values using Benjamini-Hochberg p-value correction (See also Table B6 in the appendix)
p.adjust(sapply(models, function(x) coef(summary(x))["Treatment","Pr(>|t|)"]), method = "BH")


## ------------------------------
## Figure 4
## ------------------------------

# Creating binary variables for censorship support, overall satisfaction, and willingness to protest
Experiment <- Experiment %>% 
  mutate(across(
    c(Censor_Support, Regime_Overall, Regime_Protest),
    list(Binary = ~ case_when(
      . %in% c(4,5) ~ 1,
      . %in% c(1,2,3) ~ 0,
      TRUE ~ NA_real_
    ))
  ))

# Create a data frame for the Figure 4 (group mean and standard deviations for binary measures)
Figure.4 <-  rbind(Experiment %>%
                         drop_na(Treatment) %>%
                         group_by(Treatment) %>%
                         get_summary_stats(Censor_Support_Binary, type = "mean_sd"),
                       Experiment %>%
                         drop_na(Treatment) %>%
                         group_by(Treatment) %>%
                         get_summary_stats(Regime_Overall_Binary, type = "mean_sd")
                       ,Experiment %>%
                         drop_na(Treatment) %>%
                         group_by(Treatment) %>%
                         get_summary_stats(Regime_Protest_Binary, type = "mean_sd"))
# Renaming the variables & calculating the standard errors for the group means
Figure.4 <- Figure.4 %>%
  mutate(Treatment = ifelse(Treatment == 0, "Control", "Treatment"),
         variable = case_when(
           variable == "Censor_Support_Binary" ~ "Support for Censorship",
           variable == "Regime_Overall_Binary" ~ "Regime Support: Overall",
           variable == "Regime_Protest_Binary" ~ "Willingness to Protest"
         ),
         se = sd / sqrt(n),
         ymin = mean - 1.96 * se,
         ymax = mean + 1.96 * se)
Figure.4 <- Figure.4[, c("Treatment", "variable", "mean", "ymin", "ymax")]
# Reordering the variables
Figure.4$variable <- factor(Figure.4$variable,
                                levels=c("Support for Censorship", "Regime Support: Overall",
                                         "Willingness to Protest"),
                                ordered=TRUE)
# Adjust positions between bars
pd = position_dodge2(width = 0.3, reverse = TRUE)
# left panel: censorship support
Figure.4.left <- ggplot(Figure.4[1:2,], aes(x=Treatment, y = mean, shape = variable, lty = variable)) + 
  # Point estimates
  geom_point(position = pd, size = 3, show.legend = FALSE) + 
  # Confidence intervals
  geom_linerange(aes(x = Treatment, 
                     ymin = ymin,
                     ymax = ymax),
                 size = 0.5,
                 position = pd, show.legend = FALSE) +
  # Left panel title
  ggtitle("Censorship Apparatus") +
  # Labels
  labs(y="Pr(Support)", x="") +  
  # Nicer theme 
  theme_bw() +
  scale_y_continuous(limits = c(0.45,0.65)) +
  theme(axis.text = element_text(size = 12, face = "bold"),
        axis.title = element_text(size=12),
        axis.text.y = element_text(size = 12)) 
# Middle panel: censorship support
Figure.4.middle <- ggplot(Figure.4[3:4,], aes(x=Treatment, y = mean, shape = variable, lty = variable)) + 
  # Point estimates
  geom_point(position = pd, size = 3, show.legend = FALSE) + 
  # Confidence intervals
  geom_linerange(aes(x = Treatment, 
                     ymin = ymin,
                     ymax = ymax),
                 size = 0.5,
                 position = pd, show.legend = FALSE) +
  # Middle panel title
  ggtitle("Regime: Overall Satisfaction") +
  # Labels
  labs(y="Pr(Support)", x="") +  
  # Nicer theme 
  theme_bw() +
  scale_y_continuous(limits = c(0.6,0.8)) +
  theme(axis.text = element_text(size = 12, face = "bold"),
        axis.title = element_text(size=12),
        axis.text.y = element_text(size = 12)) 
# Right panel: censorship support
Figure.4.right <- ggplot(Figure.4[5:6,], aes(x=Treatment, y = mean, shape = variable, lty = variable)) + 
  # Point estimates
  geom_point(position = pd, size = 3, show.legend = FALSE) + 
  # Confidence intervals
  geom_linerange(aes(x = Treatment, 
                     ymin = ymin,
                     ymax = ymax),
                 size = 0.5,
                 position = pd, show.legend = FALSE) +
  # Right panel title
  ggtitle("Participate in Protest") +
  # Labels
  labs(y="Pr(Willing)", x="") +  
  # Nicer theme 
  theme_bw() +
  scale_y_continuous(limits = c(0.3,0.5)) +
  theme(axis.text = element_text(size = 12, face = "bold"),
        axis.title = element_text(size=12),
        axis.text.y = element_text(size = 12)) 
# Combining the three panels
Figure.4.combined <- grid.arrange(Figure.4.left, Figure.4.middle, Figure.4.right, ncol=3)
# Saving the ggplot
ggsave(Figure.4.combined,
       file = "Figures_Paper/Figure_04.pdf",   # The directory you want to save the file in
       width = 10, # The width of the plot in inches
       height = 3.5) # The height of the plot in inches


## ------------------------------
## Figure 5
## ------------------------------

# Create a data frame for the Figure 5 (group mean and standard deviations for all three groups in study 2)
Figure.5 <-  rbind(Study2 %>%
                         drop_na(Group) %>%
                         group_by(Group) %>%
                         get_summary_stats(Censor_Support, type = "mean_sd"),
                   Study2 %>%
                         drop_na(Group) %>%
                         group_by(Group) %>%
                         get_summary_stats(Regime_Overall, type = "mean_sd"),
                   Study2 %>%
                         drop_na(Group) %>%
                         group_by(Group) %>%
                         get_summary_stats(Regime_Protest, type = "mean_sd"))
# Renaming the variables & calculating the standard errors for the group means
Figure.5 <- Figure.5 %>%
  mutate(Treatment = case_when(
    Group == "Blank" ~ "N",
    Group == "Control" ~ "P",
    Group == "Treatment" ~ "P & NP"
  ),
  variable = case_when(
    variable == "Censor_Support" ~ "Support for Censorship",
    variable == "Regime_Overall" ~ "Regime Support: Overall",
    variable == "Regime_Protest" ~ "Willingness to Protest"
  ),
  se = sd / sqrt(n),
  ymin = mean - 1.96 * se,
  ymax = mean + 1.96 * se)
Figure.5 <- Figure.5[, c("Treatment", "variable", "mean", "ymin", "ymax")]
# Reordering the variables
Figure.5$variable <- factor(Figure.5$variable,
                                levels=c("Support for Censorship", "Regime Support: Overall",
                                         "Willingness to Protest"),
                                ordered=TRUE)
# Adjust positions between bars
pd = position_dodge2(width = 0.3, reverse = TRUE)
# Left panel
Figure.5.left <- ggplot(Figure.5[1:3,], aes(x=Treatment, y = mean, shape = variable, lty = variable)) + 
  # Point estimates
  geom_point(position = pd, size = 3, show.legend = FALSE) + 
  # Confidence intervals
  geom_linerange(aes(x = Treatment, 
                     ymin = ymin,
                     ymax = ymax),
                 size = 0.5,
                 position = pd, show.legend = FALSE) +
  # Left panel title
  ggtitle("Support: Censorship Apparatus") +
  # Labels
  labs(y="Group Mean", x="") +  
  # Nicer theme 
  theme_bw() +
  scale_y_continuous(limits = c(3.3,3.7)) +
  theme(axis.text = element_text(size = 12, face = "bold"),
        axis.title = element_text(size=12),
        axis.text.y = element_text(size = 12)) 
# Middle panel
Figure.5.middle <- ggplot(Figure.5[4:6,], aes(x=Treatment, y = mean, shape = variable, lty = variable)) + 
  # Point estimates
  geom_point(position = pd, size = 3, show.legend = FALSE) + 
  # Confidence intervals
  geom_linerange(aes(x = Treatment, 
                     ymin = ymin,
                     ymax = ymax),
                 size = 0.5,
                 position = pd, show.legend = FALSE) +
  # Middle panel title
  ggtitle("Regime: Overall Satisfaction") +
  # Labels
  labs(y="Group Mean", x="") +  
  # Nicer theme 
  theme_bw() +
  scale_y_continuous(limits = c(3.6,4.0)) +
  theme(axis.text = element_text(size = 12, face = "bold"),
        axis.title = element_text(size=12),
        axis.text.y = element_text(size = 12)) 
# Right panel
Figure.5.right <- ggplot(Figure.5[7:9,], aes(x=Treatment, y = mean, shape = variable, lty = variable)) + 
  # Point estimates
  geom_point(position = pd, size = 3, show.legend = FALSE) + 
  # Confidence intervals
  geom_linerange(aes(x = Treatment, 
                     ymin = ymin,
                     ymax = ymax),
                 size = 0.5,
                 position = pd, show.legend = FALSE) +
  # Right panel title
  ggtitle("Willingness to Protest") +
  # Labels
  labs(y="Group Mean", x="") +  
  # Nicer theme 
  theme_bw() +
  scale_y_continuous(limits = c(2.9,3.3)) +
  theme(axis.text = element_text(size = 12, face = "bold"),
        axis.title = element_text(size=12),
        axis.text.y = element_text(size = 12)) 
# Combining the three panels
Figure.5.combined <- grid.arrange(Figure.5.left, Figure.5.middle, Figure.5.right, ncol=3)
# Saving the ggplot
ggsave(Figure.5.combined,
       file = "Figures_Paper/Figure_05.pdf",   # The directory you want to save the file in
       width = 10, # The width of the plot in inches
       height = 3.5) # The height of the plot in inches


## ------------------------------
## Alternative Explanations
## ------------------------------

# In the "Alternative Explanations and Limitations" section, I claim that:

# To obtain a better sense of whether the normalization theory explains the treatment
# effects above, in study 1, I asked respondents an additional question about whether they
# think “government control of the Internet is normal.” I find individuals in the treatment
# group are significantly more likely to believe that censorship is normal (β = 0.243, p = 0.002),
# and this belief is a significant mediator between the treatment and support for censorship
# (ACME = 0.149, p = 0.002), directly supporting my normalization theory.

# The numbers in this paragraph come form the following analysis: 

# Difference in belief that censorship is normal
alternative.1 <- lm(Censor_Normal ~ Treatment, data = Experiment)
summary(alternative.1)
# Causal Mediation Analysis
mediation.data.1 <- na.omit(Experiment[,c("Censor_Support","Censor_Normal","Treatment")])
mediation.1.m <- lm(Censor_Normal ~ Treatment, data = mediation.data.1)
mediation.1.y <- lm(Censor_Support ~ Treatment + Censor_Normal, data = mediation.data.1)
mediation.1 <- mediate(mediation.1.m, mediation.1.y, sims = 1000, treat = "Treatment", mediator = "Censor_Normal")
summary(mediation.1)

# To further distinguish the normalization theory from the bundling argument, in study 2,
# I included two additional outcomes measuring support for the censorship of political content
# and non-political content respectively. First, comparing the blank control and control groups,
# I find that respondents in the control group with exposure to political censorship express
# significantly lower support for not only the censorship of political content (β = −0.126,
# p = 0.004) but also non-political content (β = −0.095, p = 0.035). Then, comparing the
# control and treatment groups, I find respondents in the treatment group express significantly
# higher support for the censorship of non-political content (β = 0.101, p = 0.026). However,
# albeit in the correct direction, the effect on political content (β = 0.064, p = 0.155) is
# not statistically significant at the conventional level. Moreover, I asked whether they think
# censorship targets political content, a proxy for the perceived level of repressiveness, and also
# find it to be a significant mediator both between the blank control and the control (ACME
# = −0.017, p = 0.004), as well as between the control and the treatment at the 0.1 level
# (ACME = 0.011, p = 0.084).

# The numbers in this paragraph come form the following analysis:

# Difference in support for censorship of political/non-political
# Control vs. Blank Control
alternative.2 <- lm(Censor_Political ~ Group, data = Experiment[Experiment$Group %in% c("Blank", "Control"),])
summary(alternative.2)
alternative.3 <- lm(Censor_NonPolitical ~ Group, data = Experiment[Experiment$Group %in% c("Blank", "Control"),])
summary(alternative.3)
# Treatment vs. Control
alternative.4 <- lm(Censor_Political ~ Treatment, data = Experiment)
summary(alternative.4)
alternative.5 <- lm(Censor_NonPolitical ~ Treatment, data = Experiment)
summary(alternative.5)
# Causal Mediation Analysis: whether censorship targets political content
# Control vs. Blank Control
mediation.data.2 <- na.omit(Experiment[Experiment$Group %in% c("Blank", "Control"),c("Censor_Support","Censor_Content","Group")])
mediation.2.m <- lm(Censor_Content ~ Group, data = mediation.data.2)
mediation.2.y <- lm(Censor_Support ~ Group + Censor_Content, data = mediation.data.2)
mediation.2 <- mediate(mediation.2.m, mediation.2.y, sims = 1000, treat = "Group", mediator = "Censor_Content")
summary(mediation.2)
# Treatment vs. Control
mediation.data.3 <- na.omit(Experiment[,c("Censor_Support","Censor_Content","Treatment")])
mediation.3.m <- lm(Censor_Content ~ Treatment, data = mediation.data.3)
mediation.3.y <- lm(Censor_Support ~ Treatment + Censor_Content, data = mediation.data.3)
mediation.3 <- mediate(mediation.3.m, mediation.3.y, sims = 1000, treat = "Treatment", mediator = "Censor_Content")
summary(mediation.3)
